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Abstract 

Nonequilibrium Green's functions represent a promising tool for describing central nuclear reactions. Even at the single- 
particle level, though, the Green's functions contain more information that computers may handle in the foreseeable future. 
In this study, we investigate whether all the information contained in the Green's functions is necessarily relevant when 
describing the time evolution of nuclear reactions. For this, we carry out mean-field calculations of slab collisions in one 
dimension. 
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1. Simulations of Nuclear Reactions 

Historically, the description of central nuclear re- 
actions has involved a handful of methods. On the 
one hand, the time-dependent Hartree-Fock (TDHF) 
approach has been used to study low-energy reac- 
tions [1,2,3]. On the other hand, the semiclassical 
Boltzmann equation (BE) has been applied to the 
reactions at intermediate and high-energies [4]. Molec- 
ular dynamics approaches, sharing elements of both 
TDHF and BE, have also been employed [5,6]. Ei- 
ther TDHF or BE have some serious limitations that 
are principally remedied within the nonequilibrium 
Green's function method [7,8], as will be indicated in 
the following. Even at the basic single-particle level, 
however, the Green's functions method involves the 
handling of a vast amount of information, likely to 
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overwhelm the capabilities of computing systems. 
Here, we investigate whether all the information in the 
Green's functions is equally important for the reaction 
dynamics. To this end, we study the dynamics in one 
dimension in the absence of correlations. These results 
are principally equivalent to one-dimensional TDHF 
calculations, but are performed in terms of Green's 
functions rather than single-particle wavefunctions. 

Within the TDHF method, the system wavefunc- 
tion <1> is approximated in terms of one Slater determi- 
nant of single-nucleon wavefunctions, 

o- fc = l 

(1) 

The single-nucleon wavefunctions follow wave equa- 
tions in terms of a self-consistent mean field U: 

= + (2) 

The first applications of TDHF in the nuclear field 
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Fig. 1. Contour plots for nucleon density in a head-on col- 
lision of i^O -I- 22 Ne at Ecm = 95 MeV, from TDHF calcu- 
lations of Ref . [3] . 

have included fusion reactions. An example of the head- 
on reaction of ^''O + ^^Ne, at Ecm = 95 MeV, from 
a calculation by Umar and Oberacker [3], is shown 
in Fig. 1. In this particular collision, the nuclei form 
a compound object that subsequently breaks up, i.e. 
no fusion occurs. TDHF calculations, in fact, predict 
a so-called low-^ window, regarding fusion reactions, 
that develops at high collision energies [2,9], see Fig. 2. 
While nuclei fuse in the more peripheral reactions (at 
high angular momenta), they fail to do so in the more 
central reactions (at low angular momenta). No evi- 
dence of such a low-£ window has been found experi- 
mentally [10]. The failure of TDHF with that respect 
has been attributed to the lack of dissipation, associ- 
ated to the missing correlations in a mean-field descrip- 
tion. In this particular context, the interest in TDHF 
as a general theoretical method for addressing central 
reactions, has seriously waned. 

Higher-energy central nuclear reactions have been 
commonly described [4,11] in terms of the BE for the 
evolution of the phase-space distribution functions, 
f{r,p,t), of nucleons and other particles: 

dt dp Or dp dp ^ ' 

Here, ijjp{r,t) is the quasiparticle energy for a nu- 
cleon with momentum p, at location r, and / is a 
collision integral. BE may be solved following test- 
particle method, where the phase-space distribution 
is represented in terms of test-particles at phase-space 
locations {ri{t),p-{t)), 

f{r,p, t) ~ W ^ 5(r - r.,{t)) 5(p - pM) ■ (4) 
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Fig. 2. Fusion region in the plane of angular momentum 
vs center-of-mass energy for a -"^^O + ^^Mg reaction, from 
the TDHF calculations of Ref. [9]. The outer boundary of 
the region is associated with the fact that nuclei, in their 
relative motion, must overcome combined Coulomb and 
centrifugal barriers. 

The test-particle locations obey Hamilton-type equa- 
tions that follow from integrating the l.h.s. of Eq. (3), 

"^'^ dp ' ^'^ dr ■ 

Moreover, the test-particles undergo collisions which 
accomplish a Monte-Carlo integration of the collision 
integral / [4,11]. 

BE has been fairly successful in describing many as- 
pects of higher-energy reactions, see e.g. Fig. 3. How- 
ever, the use of BE in reactions has been criticized 
on theoretical grounds. BE relies on the quasiparticle 
picture and simple estimates [8] indicate that particle 
scattering rates are comparable to particle energies, 
which undermines that picture. In addition, in this con- 
text it is theoretically difficult [12] to separate coUi- 
sional effects, described with cross-sections and enter- 
ing the integral I, from mean-field effects entering the 
quasiparticle energies, such as in cop — /2m + U. 



2. Kadanoff-Baym Equations 

With quantum nonequilibrium many-body theory, 
the dynamics of a system, starting from an initial 
state 1$), may be described, in a self-contained man- 
ner, in terms of a generalized single-particle Green's 
function 
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Fig. 3. Proton spectra from an 800 McV/nucleon ^-^C + -'^C 
reaction. Dots represent data of Ref. [13] and histograms 
represent BE calculations of Ref. [11]. 

iG(l,l') = (*|T{V'(l)V^(l')}l$)- (6) 

Here, T is a generalized time-ordering operator which 

allows for either ordering of the singlo-particlc oper- 
ators ip and ipK The arguments represent space and 
time, 1 = {ri,ti) and discrete indices are ignored in 
the following. The generalized Green's function satis- 
fies an integral Dyson equation, 

G = Go + Go^G, (7) 

in terms of the noninteracting Green's function Go and 
the self energy E. The latter is given by 

iE(l,l') = (*|T{j(l)/(l')}l*>i.red, (8) 

where j is the source for the field tp, 

Application of the operator G^^ , inverse to Go 
in space-time, to both sides of the Dyson equation, 
yields an integro-difiierential equation of motion for 
the Green's function G. For the specific order of op- 
erators in the expectation value for G, such as in 
— iG^(l, 1') = (V'^(l') that equation reduces to 

the Kadanoff-Baym (KB) equations, 

= J dl" E+(l, 1") G*(l", 1') (10) 
+ [ dl"ES(l,l")G-(l",l'). 



Prom the two Green's functions, G^, —iG'~ yields the 
1-particle density matrix when computed at ti — ti/, 
i.e. it can be used to calculate all the 1-body observ- 
ables. For interactions of exclusively 2-body type, the 
function — iG^ can be further used to calculate the net 
system energy. 

A variety of physics phenomena, within different sit- 
uations, can be described by employing consistent ap- 
proximations to E and the dynamics. When Re E^ ~ 
U > E>, the mean-field (TDHF) approximation ap- 
plies, for which the Green's function may be approxi- 
mated by 

-iG<(l,l')«X^ 0,(1)0* (I')- (11) 

When scales associated with space-time variation in 
the average of the Green's function arguments, (1 -|- 
l')/2, are much larger than the scales associated with 
variation in the difference of the arguments, (1 — 1'), 
the quasiparticle approximation applies. Within that 
approximation, the Green's function G^ may be rep- 
resented as 

-iG<(l,l')« J dp/(p, i)e*P(='i-»'i')-*'^p(*i^*i') , 

(12) 

and a substitution of such a representation to the KB 
equations yields the BE. 

Within nuclear physics, the nonequilibrium Green's 
functions and KB equations have been primarily used 
for derivations [11,14,15,16]. In one early case, KB 
equations have been solved directly for a uniform sys- 
tem of equilibrating interpenetrating nuclear matters 
[8], see Fig. 4. The system of nuclear matters repre- 
sents the first stages of a central collision of nuclei 
at intermediate incident energies. The KB results for 
this model uniform system have been compared to 
BE calculations [8,17]. Differences in these results, 
see Fig. 4, in the way in which the momentum space 
gets populated, and in the pace of approach to equi- 
librium, demonstrate the importance of going beyond 
the quasiparticle approximation in the description of 
intermediate-energy central reactions. More results 
from solving directly the KB equations have been 
obtained in other fields [18,19,20]. 



3 




Fig. 4. Contour plots of the distribution of nucleons in mo- 
mentum, for the evolution of interpenetrating nuclear mat- 
ters that start, at t = 0, from two separated Fermi spheres 

with displacement corresponding to the incident energy of 
^lab = 400 McV/nuclcon. The left column displays the evo- 
lution for BE and the center and right columns display 
the evolution for KB equations when starting, respectively, 
from an uncorrelated and correlated initial states, after 
Ref. [8]. 

3. Towards Direct Application of the Green's 
Functions to Reactions 

The Green's function technique is attractive in that 
it allows one to incorporate different types of correla- 
tions within E in the system dynamics, without resort- 
ing to the quasiparticle limit. The description of the 
reaction can be made consistent with the description 
of the structure of the reacting nuclei [21,22]. The se- 
rious disadvantage of the method is that potentially 
enormous amount of information, representing the 8- 
dimensional structure of the Green's functions, must 
be monitored for a reaction. A secondary deficiency 
is that potentially two separate calculations need to 
be carried out to complete a reaction simulation: one 
preparing the system for the reaction and second ad- 
vancing the reaction in time. 

The 8-dimensional information in the Green's func- 



tion can be generally understood as mappable onto 

the 3 components of particle momentum, onto energy 
(an independent variable outside of the quasiparticle 
limit) and onto time and 3 components of position 
where the particle momentum and energy arc consid- 
ered. An increase in the dimension by one variable, the 
energy, as compared to BE, might be considered a mi- 
nor complication. However, the serious issue is that no 
positive-definite quantity can be associated with the 
Green's functions, such as /, thus precluding the pos- 
sibility of applying the test-particle method. The ad- 
vantage of the test-particle method is that it allows 
one to concentrate on the parts of the phase space that 
matter at a given stage of the reaction. Without this 
method, a mesh in space-time needs to be established 
and diflferent space-time regions may potentially need 
to be monitored, irrespectively of whether or not those 
regions dominate a given stage of a reaction. An alter- 
native to the mesh is projection onto mean-field wave- 
functions [23], but that can limit excitation energies 
in the reaction simulations. In the following, we shall 
try to assess whether the whole involved space-time 
region needs indeed to be monitored within the mesh 
set up for the reaction. We shall carry out the assess- 
ment in one dimension, where the cost of monitoring 
the whole space-time is moderate, unlike in three di- 
mensions, where the cost of such monitoring is devas- 
tating. We shall also consider the possibility of prepar- 
ing the initial state when following time evolution such 
as in a reaction. 

3.1. Slabs in One Dimension 

To assess the importance of different regions in the 
domain of the Green's functions, we consider interac- 
tions of nuclear slabs, for which density changes in the 
x-direction and for which transverse degrees of freedom 
are frozen. Those slabs arc spin-isospin symmetric and 
follow mean-field dynamics with Coulomb interactions 
switched off. Under these circumstances, the dynamics 
is consistently described in terms of the Green's func- 
tion at equal time arguments and at equal trans- 
verse positions. The trace of the density matrix over 
spin and isospin is then 

p{x,x',t) = -AG^{x,y,z,t,x' ,y,z,t), (13) 

where the factor of 4 represents the spin-isospin de- 
generacy. The density matrix satisfies the equation of 
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motion 



i-Q-^p{x,x' ,t) = 



p{x,x',t) 



2m dx'' 



+ Uix',t) 



p{x,x',t) , 
(14) 



where U is & Skyrme mean-field of the form 



3 , , 2 + 0" 
(7 = - to n[x, t) H — t3 [n(x, t)\ 



(15) 



and the density n represents diagonal elements of the 
density matrix, n{x, t) = p(x, x, t). The parameter con- 
stants in the Skyrme interaction are fitted to the prop- 
erties of nuclei. If the density matrix is diagonalized 
in terms of eigenstate wavefunctions, as in (11), then 
Eq. (14) can be solved by having the wavefunctions fol- 
low the TDHF equations (2). However, in the context 
of the Green's functions, it is the equation for the den- 
sity matrix, Eq. (14), that is of interest. Net kinetic 
energy associated with the x-direction, per unit trans- 
verse area, is 



Ax ■ 



— n(x, x\t)\ , 



2m J dxdx' 
and the net potential energy per unit area is 



V{t) 



dx 



-to [n{x)] + —t3 [n(x)] 



2+cr 



16 



(16) 



(17) 



In evolving the slab systems, we use periodic bound- 
ary conditions for the calculational region and employ 
the split-operator method to advance the density ma- 
trix in time. Thus, the change in the density operator 
over the time-interval At may be represented as 



p{x,x' , t + At) 



(18) 



where /C is the kinetic-energy operator. Using the 
Baker-Campbell-Hausdorff formula, the exponentials 
of a sum of non-commuting operators in (18) can be 
rewritten into a product of exponentials of individual 
operators 



J{K+U) At 



e^^'^e^'^^'e^'^'r +o[{Atf] 



(19) 



with the accuracy of the second order in At. The split- 
operator method refies on the above rearrangement 
and the fact that the operators K. and U axe diago- 
nal either in the momentum or configuration space. 
Application of an exponential in either space amounts 
then just to a multiplication by a phase factor. Within 
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Fig. 5. Evolution of the energy per nucleon (top) and of 
the size of density distribution (bottom) when adiabati- 
cally switching the potential from a harmonic oscillator to 
a mean-field one, for different values of the r parameter in 
Eqs. (20) and (21). 

the method, the advancement of the density matrix 
in time can be then accomplished by transforming the 
matrix between configuration and momentum spacing, 
exploiting Fast-Fourier- Transformation, and multipli- 
cation by the respective phase factors. Unitarity is built 
into the method. 

With regard to the evolution, we first consider the 
possibility of preparing the initial state for a reaction 
in an adiabatic fashion [24]. We start out with a state 
obtained by filling out two lowest orbits of the har- 
monic oscillator potential and we gradually change the 
potential from that of the harmonic oscillator to the 
mean-field one, with 

Uf{x, t) = f{t) ^x^ + [1 - fit)] U{n{x, t)) , (20) 



where we typically use 



-to)/t 



(21) 



Figure 5 shows the evolution of the net energy per nu- 
cleon and net size of the density distribution, when 
gradually changing the potential. Notably, the energy 
itself is not that good indicator of the arrival at the 
ground state, because it depends quadratically on the 
deviation from the ground state. Characteristics of the 
density, on the other hand, depend linearly on the devi- 
ation and are better indicators. Values of r > 20 fm/c 
lead to a very good approximation of the mean-field 
ground-state. 
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Fig. 6. Evolution of density for slabs colliding at cm. energy 
of 0.1 MeV/nucleon. 

We next consider slab collisions. For the collisions, 
the slabs are boosted and positioned against each 
other. In a boost, the individual density matrices 
transform according to 



p{x, x' ,t = 0) 



p{x, x' ,t = 0) e 



(22) 



where pcm is the initial cm. momentum per nucleon. 
Figure 6 shows the evolution of a slab reaction at inci- 
dent cm. energy of 0.1 MeV/nucleon. In the absence of 
Coulomb interactions, fusion occurs at this energy. At 
energies 15 MeV/nucleon > Ecm > 0.6 MeV/nucleon, 
the slabs pass through each other, but get excited 
and may break up. At still higher energies, Ecm ^ 
15 MeV/nucleon, the identity of the original slabs is 
largely lost, see Fig. 7, and the system breaks into a 
multitude of fragments. In the nuclear reaction termi- 
nology, this is called multifragmentation. 

Figure 8 shows next the intensity plot for the real 
part of the density matrix for the system of slabs col- 
liding at the cm. energy of 25 MeV/nucleon. The den- 
sity matrix is always largest along the x — x' diago- 
nal and the real components tend to marginally domi- 
nate over the imaginary components. The values of the 
matrix along the diagonal represent the local density, 
i.e. Fig. 7. Early on in the reaction, the more significant 
values of the density matrix are limited to square-like 
regions. The diagonals of these squares are coUinear 
with the diagonal of the density matrix and their size 
represents the support of the nucleon wavefunctions 
within the individual slabs. When the system fuses, 
the size of this region expands to one representing the 
size of the fused system. The system evolves further 
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Fig. 7. Evolution of density for slabs colliding at cm. energy 
of 25 MeV/nucleon. 
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Fig. 8. Intensity plots for the real part of the den- 
sity matrix p{x,x',t) in the collision of slabs at 
Ecm = 25 MeV/nucleon. The snapshots are taken at the 
same times as in Fig. 7 

and breaks up into pieces, and the region of significant 
values for the density matrix expands. Besides regions 
adjacent to the diagonal, islands and streaks of signif- 
icant values emerge, moving away from the diagonal 
with time. These are associated with the fragmenta- 
tion of the original single-nucleon wavefunctions, into 
pieces taking off with different fragments. The regions 
of significant off-diagonal values for p(x,x' ,t) repre- 
sent phase correlations between the amplitudes for the 
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Fig. 9. Manipulation of ofT-diagonal elements of the den- 
sity matrix p{x,x',t) in testing. White regions represent 
elements that get suppressed. 



nucleons, from the individual original states, to move 
away with one or another fragment. 

The need to monitor far off-diagonal values of the 
Green's functions would have been devastating for the 
capability of carrying reaction simulations in 3 dimen- 
sions. However, if the fragments from the break-up are 
not likely to ever meet again, the phase information in 
the far off-diagonal values should not matter. To test 
the importance of these far off-diagonal elements, we 
repeat the calculations employing now a strong imagi- 
nary superoperator potential that suppresses elements 
away from the diagonal axis and leaves the vicinity of 
the axis in the density matrix intact, see Fig. 9. Note 
that the periodic boundary conditions for the system 
imply a tile-like periodicity of the density matrix p 
in the x-x' plane, with the values on the diagonal re- 
peated next on the lines passing through the corners at 
(— L, L) and (L, —L). When the periodicity is imposed 
onto the superoperator, the suppression of the matrix 
elements has to occur in valleys in the x-x' plane that 
are separated by ridges, where the elements remain in- 
tact. 

Figure 10 compares the results for the density at 
three characteristic times in the slab collision at Ecm = 
25 MeV/nucleon. Each line corresponds to evolutions 
with different retained regions off the diagonal. Only 
when the matrix elements are suppressed from \x — 
x'\ = 5 fm on, some changes in the density, compared to 
the standard evolution, begin to emerge at late times. 
Even these changes are subtle and the system may still 
reach about the same final state, just slightly later. Less 
severe trimmings of the elements leave no visible signs 
in the evolution of density for times relevant for reac- 
tion dynamics. Figure 11 shows in addition the evolu- 
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Fig. 10. Density at different stages of the evolution at the in- 
cident cm. energy of 25 MeV/nucleon. Bottom curves rep- 
resent the density from the standard evolution. The other 
curves, from bottom up, show the density obtained from 
an evolution where the matrix elements at \x — x'\ > 20, 
15, 10, and 5 fm, respectively, are suppressed. For clarity, 
the results for the density of the different evolutions are 
staggered by 0.1 fm~^. 
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Fig. 11. Net energy components (top) and system size 
(bottom) for elements of density matrix suppressed at 
\x ~ x'\ > 2x0, at different values of xq. 

tion of the system size and of the net energy compo- 
nents for different off-diagonal suppressions. Only for 
elements eliminated at — > 5 fm, a minute change 
in the energy breakdown begins to emerge at late times. 

Reasoning that relative phases of amplitudes asso- 
ciated with different fragments would matter only if 
the fragments could meet again, we have further car- 
ried out calculations where we allowed the fragments 
to reach the boundary of the calculational region, cross 
to an adjacent region and collide there again. Even 
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in those cases, we found no perceptible difference in 

the evolution of the density when elements had been 
eliminated at distances |x — x'| > 10 fm. The only sit- 
uation when we find some genuine impact of the far 
off-diagonal elements is when trying to run the evolu- 
tion backwards in time. The original mean-field evo- 
lution is time reversible and such, to a very good ex- 
tent, is its numerical realization. However, reversibility 
gets lost when erasing off-diagonal elements and con- 
sidering longer time intervals. Thus, e.g. if we suppress 
elements at \x — x'\ > 10 fm, while evolving the sys- 
tem forward in time, past the reaction, and then evolve 
the system backwards with element suppression, we 
find that the system has a difficulty reaching its initial 
state [25] . In general, the system evolved backwards in 
time, with elements suppressed, acquires, in the past, 
some features of the state in the future. Thus, in the en- 
ergy region for fusion, the system is likely to stay fused 
when evolved back to the early times. In the energy 
range for multifragmentation, the system may evolve 
back to two residues but, between them, an additional 
fragment may remain and the two residues may stay 
excited. 



4. Summary 

To sum up, we have obtained encouraging results 
in our study of slab collisions in one dimension. Only 
a limited range of arguments in the density matrix is of 
practical importance for the evolution forward in time. 
The scale for that range turns out in practice to be set 
by the inverse of Fermi momentum. Further, we find 
that an adiabatic transformation of the Hamiltonian 
allows for an arrival at the ground state of an inter- 
acting, albeit mean-field, Hamiltonian. Those results 
bode well for the possibility of simulating central nu- 
clear reactions, in 3 dimensions, following the method 
of nonequilibrium Green's functions. 
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